home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / spline.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  153 lines
  1. ; $Id: spline.pro,v 1.5 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1983-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. function spline,x,y,t,sigma
  7. ;+
  8. ; NAME:
  9. ;    SPLINE
  10. ;
  11. ; PURPOSE:
  12. ;    This function performs cubic spline interpolation.
  13. ;
  14. ; CATEGORY:
  15. ;    Interpolation - E1.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;    Result = SPLINE(X, Y, T [, Sigma])
  19. ;
  20. ; INPUTS:
  21. ;    X:    The abcissa vector. Values MUST be monotonically increasing.
  22. ;
  23. ;    Y:    The vector of ordinate values corresponding to X.
  24. ;
  25. ;    T:    The vector of abcissae values for which the ordinate is 
  26. ;        desired. The values of T MUST be monotonically increasing.
  27. ;
  28. ; OPTIONAL INPUT PARAMETERS:
  29. ;    Sigma:    The amount of "tension" that is applied to the curve. The 
  30. ;        default value is 1.0. If sigma is close to 0, (e.g., .01),
  31. ;        then effectively there is a cubic spline fit. If sigma
  32. ;        is large, (e.g., greater than 10), then the fit will be like
  33. ;        a polynomial interpolation.
  34. ;
  35. ; OUTPUTS:
  36. ;    SPLINE returns a vector of interpolated ordinates.
  37. ;    Result(i) = value of the function at T(i).
  38. ;
  39. ; RESTRICTIONS:
  40. ;    Abcissa values must be monotonically increasing.
  41. ;
  42. ; EXAMPLE:
  43. ;    The commands below show a typical use of SPLINE:
  44. ;
  45. ;        X = [2.,3.,4.]      ;X values of original function
  46. ;        Y = (X-3)^2         ;Make a quadratic
  47. ;        T = FINDGEN(20)/10.+2     ;Values for interpolated points.
  48. ;                    ;twenty values from 2 to 3.9.
  49. ;        Z = SPLINE(X,Y,T)     ;Do the interpolation.
  50. ;
  51. ;
  52. ;
  53. ; MODIFICATION HISTORY:
  54. ;    Author:    Walter W. Jones, Naval Research Laboratory, Sept 26, 1976.
  55. ;    Reviewer: Sidney Prahl, Texas Instruments.
  56. ;    Adapted for IDL: DMS, Research Systems, March, 1983.
  57. ;
  58. ;-
  59. ;
  60. on_error,2                      ;Return to caller if an error occurs
  61. if n_params(0) lt 4 then sigma = 1.0 else sigma = sigma > .001    ;in range?
  62. n = n_elements(x) < n_elements(y)
  63. ;
  64. if n le 1 then message, 'X and Y must be arrays.'
  65.  
  66. ;
  67. xx = x * 1.            ;Make X values floating if not.
  68. yp = fltarr(n*2)        ;temp storage
  69. delx1 = xx[1]-xx[0]        ;1st incr
  70. dx1=(y[1]-y[0])/delx1
  71. ;
  72. nm1 = n-1L
  73. np1 = n+1L
  74. if (n eq 2) then begin
  75.     yp[0]=0.
  76.     yp[1]=0.
  77.    end else begin
  78.     delx2 = xx[2]-xx[1]
  79.     delx12 = xx[2]-xx[0]
  80.     c1 = -(delx12+delx1)/delx12/delx1
  81.     c2 = delx12/delx1/delx2
  82.     c3 = -delx1/delx12/delx2
  83. ;
  84.     slpp1 = c1*y[0]+c2*y[1]+c3*y[2]
  85.     deln = xx[nm1]-xx[nm1-1]
  86.     delnm1 = xx[nm1-1]-xx[nm1-2]
  87.     delnn = xx[nm1]-xx[nm1-2]
  88.     c1=(delnn+deln)/delnn/deln
  89.     c2=-delnn/deln/delnm1
  90.     c3=deln/delnn/delnm1
  91.     slppn = c3*y[nm1-2]+c2*y[nm1-1]+c1*y[nm1]
  92.    endelse
  93. ;
  94.     sigmap = sigma*nm1/(xx[nm1]-xx[0])
  95.     dels = sigmap*delx1
  96.     exps = exp(dels)
  97.     sinhs = .5d0*(exps-1./exps)
  98.     sinhin=1./(delx1*sinhs)
  99.     diag1 = sinhin*(dels*0.5d0*(exps+1./exps)-sinhs)
  100.     diagin = 1./diag1
  101.     yp[0]=diagin*(dx1-slpp1)
  102.     spdiag = sinhin*(sinhs-dels)
  103.     yp[n]=diagin*spdiag
  104. ;
  105.     if  n gt 2 then for i=1L,nm1-1 do begin
  106.         delx2 = xx[i+1]-xx[i]
  107.         dx2=(y[i+1]-y[i])/delx2
  108.         dels = sigmap*delx2
  109.         exps = exp(dels)
  110.         sinhs = .5d00 *(exps-1./exps)
  111.         sinhin=1./(delx2*sinhs)
  112.         diag2 = sinhin*(dels*(.5*(exps+1./exps))-sinhs)
  113.         diagin = 1./(diag1+diag2-spdiag*yp[n+i-1])
  114.         yp[i]=diagin*(dx2-dx1-spdiag*yp[i-1])
  115.         spdiag=sinhin*(sinhs-dels)
  116.         yp[i+n]=diagin*spdiag
  117.         dx1=dx2
  118.         diag1=diag2
  119.        endfor
  120. ;
  121.  
  122.     diagin=1./(diag1-spdiag*yp[n+nm1-1])
  123.     yp[nm1]=diagin*(slppn-dx2-spdiag*yp[nm1-1])
  124.     for i=n-2,0,-1 do yp[i]=yp[i]-yp[i+n]*yp[i+1]        
  125. ;
  126. ;
  127.     m = n_elements(t)
  128.     subs = replicate(long(nm1),m) ;subscripts
  129.     s = xx[nm1]-xx[0]
  130.     sigmap = sigma*nm1/s
  131.     j=0L
  132.     for i=1L,nm1 do $ ;find subscript where xx[subs] > t(j) > xx[subs-1]
  133.         while xx[i] gt t[j] do begin
  134.             subs[j]=i
  135.             j=j+1
  136.             if j eq m then goto,done
  137.             endwhile
  138.     
  139. done:    subs1 = subs-1
  140.     del1 = t-xx[subs1]
  141.     del2 = xx[subs]-t
  142.     dels = xx[subs]-xx[subs1]
  143.     exps1=exp(sigmap*del1)
  144.     sinhd1 = .5*(exps1-1./exps1)
  145.     exps=exp(sigmap*del2)
  146.     sinhd2=.5*(exps-1./exps)
  147.     exps = exps1*exps
  148.     sinhs=.5*(exps-1./exps)
  149.     spl=(yp[subs]*sinhd1+yp[subs1]*sinhd2)/sinhs+ $
  150.         ((y[subs]-yp[subs])*del1+(y[subs1]-yp[subs1])*del2)/dels
  151.     if m eq 1 then return,spl[0] else return,spl
  152. end
  153.